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We study fully occupied lattice systems of classical magnetic dipoles which point along random 
axes. Only dipolar interactions are considered. From tempered Monte Carlo simulations, we obtain 
numerical evidence that supports the following conclusions; in three dimensions, (a) there is an 
equilibrium spin glass phase at temperatures below Tc, where ksTc — (0.86 ± 0.07)ed and Ed is a 
nearest neighbor dipole-dipole interaction energy, (b) in the spin glass phase the overlap parameter 
is approximately given by — T /Tc, and (c) the correlation length ^ diverges at Tc with a critical 
exponent i/ = 1.5±0.5;in two dimensions ^ diverges at or near T = and = 3 ± 1. 
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I. INTRODUCTION 

Several decades after the experimental discovery of 
spin glasses,- convincing numerical evidence for an equi- 
librium phase transition between the paramagnetic and 
spin-glass phases of the random bond Ising^ii^ model in 
three dimensions is at last available. Somewhat more 
controversial evidence is also available for the Hcisen- 
berg model4ii^ No such results that we know of exist 
for systems in which dipole-dipole interactions dominate. 
This is in spite of all the interest that has arisen in 
these systems since nanosized magnetic particles have 
become experimentally available.— Randomness, one of 
the two essential ingredients for spin-glass behavior, can 
arise from spatial disorder,^ which in turn, most of- 
ten, brings about random magnetic anisotropics. One 
might naively expect that the long range nature of dipo- 
lar interactions would only strengthen the spin glass 
phase that is observed in the random (nearest neigh- 
bor) bond Ising model. However, recent results from 
computer simulations suggest that an equilibrium spin 
glass phase does not obtain in a spatially disordered 
system of magnetic dipoles which point along parallel 
axesi^ The reason for this somewhat unexpected result 
may be the nature of frustration that is peculiar to dipo- 
lar systems. In them, there is frustration whether they 
are spatially ordered or not. It is precisely because of 
this that ferro- or antiferro-magnetism prevails in well 
ordered crystalline dipolar systems depending delicately 
on lattice geometry.^ On the other hand, spin-glass like 
behavior has been observed in experiment a^°i^^'^^ and 
in simulation&i^iii^ii^ii^iii of dipolar systems with ran- 
dom anisotropics, but all this evidence comes from out of 
equilibrium phenomena, as exhibited by time dependent 
susceptibilities, nonexponential relaxation, and aging. 

We study the equilibrium behavior of systems of in- 
teracting magnetic dipoles which are oriented along ran- 
dom anisotropy axes in two and three dimensions (D). 
This random axes dipolar (RAD) model is like the old 
model of Harris, Plischke, and Zuckermanji^ except that 
we deal with dipole-dipole, rather than nearest neighbor 



(nn) interactions. Some motivation for the RAD model 
comes from the fact that anisotropy energies in nanopar- 
ticle assemblies are ofteni^ much larger than the dipole- 
dipole interaction energy between two nearest neighbors. 
As in an Ising model, spins in the RAD model can only 
point "up" or "down" along each one of their own axes, 
as is discussed in Refs. [Ullaill- Two independent ran- 
dom numbers per site are needed to determine all axes 
directions, which is the same number as for a nn ran- 
dom bond Ising model on a square lattice, though the 
interaction range is of course quite different. 

When we simulate the time evolution of the RAD 
model, we flip each spin up and down along its own axis. 
We thus make no attempt to simulate how each individ- 
ual spin overcomes large anisotropy barriers. Rather, we 
expect our simulations to mimic the collective time evo- 
lution effects that follow after single spin energy barriers 
are surmounted, as illustrated in Figs. 1 and 2 of Rcf. 
[l^ . Anyway, our main interest does not lie in the time 
dependent properties of the RAD model, but in its equi- 
librium behavior, which must clearly be the same as for a 
system of magnetic dipoles under a dominant anisotropy 
with random axes. 

A summary of our results follows. We first illus- 
trate advantages of tempered^-^ Monte Carlo (TMC) over 
Metropolis^^ Monte Carlo (MMC) simulations for the 
calculation of equilibrium behavior. This includes a com- 
parison of the time dependent magnetic susceptibility x 
(from MMC runs), which is characteristic of spin glasses, 
for the RAD model in 2D, with equilibrium results that 
follow from TMC simulations. We obtain equilibrium 
results (from TMC simulations) for systems of L*^ spins 
{d is the lattice dimension) for d = 2 and d = 3, for 
L = 4, 8, 16 and for L = 4,6, 8, 12, respectively. Simula- 
tions of larger systems are very time consuming, because 
running times grow as L^'^ for systems with dipolar in- 
teractions. Extrapolations to the L — > oo limit point 
to the following conclusions. In three dimensions (3D), 
the paramagnetic phase covers the T > Tc range, where 
T is the temperature, Tc = (0.86 ± 0.07)ed/fcs, fcs is 
Boltzmann's constant, and Sd is a dipole-dipole nn inter- 
action energy which is defined below, in Sec. Ill Al For 
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T < Tc, there is an equilibrium spin glass phase. In it, 
the overlap parameter, as defined in Sec. Ill C[ is approxi- 
mately given by ^1 — T jT^. From our results we cannot 
quite conclude whether the droplet--'.?^ qj. j^gg25,26 pjg_ 

ture describes the RAD spin glass in 3D. Results for the 
correlation length ^, are consistent with C ~ (^^ ~ 2^c) 
where ~ 0.88 and = 1.5 ± 0.5.^"^ In 2D, the param- 
agnetic phase covers the T > range, though we cannot 
rule out a spin glass phase below T ~ 0.1. Results for the 
correlation length ^ are consistent with ^ ~ T~'^, where 
iy = 3±l. 



II. MODEL AND METHOD 
A. Model 

To define the model, let 

ij Q/3 

be its Hamiltonian, where 5*^ is the a (one of three) 
component of the classical spin on a cubic lattice site i, 

Tf - e4a/n,nSc.p - ir^r^lrl), (2) 

Tij is the distance between i and j, is an energy, and a 
a nn distance. Each spin points along a randomly chosen 
direction. More precisely, let Uj be a 3— component vec- 
tor chosen randomly for each % from a spherically uniform 
distribution of unit vectors, and let Uj = ±1 at each site, 
such that Sj — UjCFj. Then, Ti. becomes, 

T-i = --^J^Ja^aj, (3) 

ij 

where Jy = - Ea,/3 '^i^^<"f ■ Thus, the RAD model 
is an Ising model whose bonds Jij are determined by 
the dipole-dipole terms T^"'' and the set of 3-component 
randomly oriented unit vectors {uj}. 

We use periodic boundary conditions in 2D and 3D. 
Simple cubic lattices and zero applied magnetic field H 
are assumed throughout. We only work with L"^ box- like 
systems, and let dipole-dipole interactions act between 
each spin and all other spins within an L"^ box centered 
on it. Because of the long range nature of dipolar interac- 
tions, contributions from beyond this box would have to 
be taken into account (by some scheme, such as Ewald's 
summation) if spins were to point in any one preferred 
direction. They do not do so in this (nonferromagnetic) 
model as long as — 0. The boundary conditions as 
well as the L'^ box scheme we use here are as in Refs. 
diilll. Finally, it is worth recalling that thermal equi- 
librium results obtained for H=0 for large cubic-shaped 
systems can, by virtue of Griffiths theorem^ be general- 
ized to other shapes in three dimensions. 

From here on, all temperatures are given in terms of 
Sd/^B, where ks is Boltzmann's constant. 



B. Monte Carlo 

Let us first specify how we update the state of the sys- 
tem in all Monte Carlo evolutions. Initially, we compute 
the dipolar field at each site. Throughout a computer 
run, tables of all spins and dipolar fields are kept. Dipo- 
lar fields are updated throughout all sites in the system 
whenever a spin is fiipped. Thus, no computer time is 
wasted whenever an attempt to flip a spin ends in fail- 
ure. This becomes important at low temperatures. 
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FIG. 1: (Color online) qf and 52 v s time for systems of 6 x 6 x 6 
spins at T = 0.5. □ and o are for 52 and ql, respectively. They 
both follow from the MMC algorithm. On the other hand, ■ 
and • are also for 52 and ql , respectively, but they both follow 
from the TMC algorithm. All data points stand for averages 
over 200 samples, each with different random anisotropy axes. 
All systems were allowed to evolve for over 10^ MCS before 
any measurements were taken. 

The idea behind the tempered Monte Carlo 
algorithmf22, is as follows. Consider two systems, 2 
and 1, in thermal equilibrium, not among themselves 
but each one of them with its own heat bath, at 
temperatures Ti and T2, respectively. Let T2 > Ti, and 
let E2 and Ei be the energies of systems 1 and 2 at 
some given time. In the TMC algorithm, the states of 
two systems are exchanged with a certain probability 
p at some specified times. It follows straightforwardly 
that the canonical thermal probability distributions for 
systems 1 and 2 are unchanged ii p = 1 ii E2 < Ei, 
and p — exp[(/3i — l32){Ei — E2)] if E-? > Ei, where 
= 1/Tfc for fc = 1,2.22 

We do TMC simulations on k identical systems at tem- 
peratures T + nAT, where n — 1, 2 . . . fc, with initially 
independent random spin configurations, and let all sys- 
tems evolve in time following the MMC algorithm for 
a number n of consecutive MMC sweeps. (We choose 
71 = 10 throughout.) After every n MMC sweeps, pairs 
of systems are given a chance to exchange energy, follow- 
ing the above given rule. More specifically, systems 2n 
and 2n + 1, for n — 0, 1, 2, • • • , are allowed to exchange 
states at jn times, where j = 1, 2, • • • , and systems 2n 
and 2n — 1, for n = 1, 2, • • • , are allowed to exchange 
states at {j + l/2)h times. We choose AT as follows. Fre- 
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FIG. 2: (Color online) (a) In plane (x||) and out of plane (x±) 
susceptibilities vs T for 2D systems oi L x L spins. A, I, and 
• are for X||, from MMC runs of 10*, 1Q^ and lO'' MCS, 
respectively. The low lying string of ♦ symbols is for x±- B 
is for data points which follow from TMC for X||- Finally, all 
ffl. A, and ♦ are for L — 32 and the rest of data points are for 
L = 16. Error bars are roughly given by icon sizes, (b) 1/x 
vs T, from TMC simulations. All data points above (below) 
the diagonal dotted line are for 3D (2D) systems, o, □, o, and 
A are for L = 4, 8, 16, and 32, respectively. ♦, EB, ■, and T 
are for L = 4, 6, 8, and 12, respectively. Parameters for TMC 
runs are given in table I. 

quent exchanges take place if the energy difference AE 
between systems 2n and 2n ± 1 is not much larger than 
the energy fluctuationsjiS^i On the other hand, we know 
from our own simulations of the RAD model, that the 
specific heat C fulfills C w for T < 0.6 and C < 
for all < T, both in 2D and 3D. Using C < T^, one 
obtains AT < 1/^/N, which is the desired condition. 

How much faster stationary states are approached in 
TMC than in MMC simulations is illustrated in Fig. [TJ 
where plots of qf and q2 (defined in Sec. Ill Cp vs time are 
shown, using data points from both MMC and TMC. For 
further comparison, results obtained for the susceptibility 
X from MMC and TMC simulations are shown in Fig. 
[5^. All data points, except the low lying branch, are 
for the "in plane" susceptibility X||> that is, the linear 
in plane magnetization response to an in plane applied 
magnetic field. The low lying branch in Fig. is for 
the "out of plane" linear susceptibility x±- (As is well 
known, dipolar interactions lead to "shape anisotropy" , 
which for 2D gives X|| ^ X±-^") We often write x for 
X||- All data points in Figs. [2^, and[2J: follow from 
measurements of magnetization fluctuations in _ff = 0. 



The data points from MMC simulations clearly ex- 
hibit time dependent effects that are sometimes associ- 
ated with spin glasses. The peak in x\\ shifts to lower 
values of T as the number of MCS increases. This is as 
expected from a spin glass. Results from MMC in 3D 
(not shown) do not differ qualitatively from the results 
shown in Fig. [2^ for 2D. Finally, equilibrium susceptibil- 
ities that follow from magnetization fluctuations in TMC 
runs are shown in Fig. [2)3 for system of various sizes, in 
2D and in 3D. 



C. Overlaps 

We next define the equilibrium quantities we calculate. 
Following the original idea of Edwards and Anderson, '^^ 
consider two identical replicas, 1 and 2, of a system. Both 
replicas have the same set of anisotropy axes but evolve 
in time independently, starting from arbitrarily different 
initial statesi^ Let 

'/'.^'^fVf, (4) 

where aj^'' and crj^'' be the spins on site j of replicas 1 
and 2, and 

q = L-''Y,cl>,. (5) 

i 

We also define the moments of q, qk = {\ q \'') , ioi: k — 1, 
2 and 4, where (. . .) stands for an average over equilib- 
rium states of a large number Ng of replica pairs with 
independent random axes orientations. Note we use an 
absolute value in the definition of qi. We refer to qi as 
the overlap parameter. The spin glass susceptibility is 
given by L'^q2- 

Recall that if the probability distribution P{q) in the 
spin glass phase differs from zero only in a vanishingly 
small neighborhood of some q = ±go, where < qo < 
as in the droplet modelfSliS^ then, q2 ^ ql > 0. On the 
other hand, if P{q) ^ over a finite range of q values, as 
in the RSB schemer^»2£ then q2 > qf. 

In order to keep track of time evolutions, we also de- 
fine (f>j(tQ,t) — <Tj{to)aj(to +t), in close analogy with the 
definition of Eq. 21 except that both aj{to) and aj{to+t) 
are the same spin, at site j, at times to and to + t, re- 
spectively. We also define q{t,to) = L~''-^ - (f)j{to,t), and 
the moments qk in obvious analogy to qk. No measure- 
ment is ever taken, neither for the calculation of qk nor 
for qk{to,to +t), in any simulation up to time to. The 
question is how to choose to- Obviously, the i — > oo 
limit of qk{t,to) depends on to- Indeed, aging is the out- 
come of a rather long lasting dependence on tp'^'^'^ 
For equilibrium results, we choose sufficiently large val- 
ues of to in order that qk{t,to) reach steady state before 
t = to- Failure to do so would imply that equilibrium 
had not been reached by to^ after which time measure- 
ments had been taken. We thus (a) let to be halfway to 
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TABLE L TMC simulation parameters. 



d, L 


2, 4 


2, 8 


2, 16 


2, 32 


3, 4 


3, 6 


3, 8 


3, 12 


Ns 


1000 


600 


300 


100 


1800 


800 


400 


175 


AT 


0.1 


0.1 


0.1 


0.05 


0.1 


0.1 


0.05 


0.05 


MCS 


10^ 


10^ 


10^ 


10^ 


4 X 10^ 


2 X 10^ 


2 X 10^ 


2 X 10^ 



the end of each MC run, that is, we let to = tf, where 
2t f is the total number of MC sweeps taken in any given 
run, starting from a random spin configuration, and (b) 
let tf he sufhcicntly large for qk to have reached steady 
state by the end of the run. For short, we write qk for 
qk{tf /2,tf). All of this is necessary but not sufficient. 
Conceivably, an exceedingly fast initial evolution away 
from a disordered state at an early to could drive gfe(to) 
to a null value, long before equilibrium was reached. On 
the other hand, the value of qk, averaged over the time 
interval (to,t/), would still depend ontf. Therefore, for 
equilibrium calculations we choose to (and therefore tf) 
sufficiently long for qk to become equal to qk- For com- 
parison, equilibrium data points for both qk and qk are 
sometimes displayed jointly. 

III. EQUILIBRIUM RESULTS 

We report our equilibrium results in this section. The 
relevant parameters for all TMC simulations from which 
these results follow can be fomid in Table I. 




FIG. 3: (Color online) Log-log plots of 52 (black) and (red) 
vs 1/L in 3D. Closed and open icons are for 52 and qf, re- 
spectively. ♦ and o are for T = 0.45, • and o are for T — 0.6, 
■ and □ are for T — 0.8, T and V are for T = 1.0, and A and 
A are for T = 1.1. Lines are guides to the eye. 

Plots of equilibrium values of q2 and qf vs 1/L are 
shown in Fig. [3] for systems in 3D at various temper- 
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FIG. 4: Plots of q2 vs 1/L in 2D for T=0.2 (♦), 0.4 (a), 0.6 
(T), 0.8 (■), 1.0 (•). Lines are guides to the eye. 

atures. For T < 1 {T > 1), q2 and qf curve upward 
(downward). This suggests Tc ^ 1. Extrapolations per- 
formed on linear plots (not shown) of 52 and qf vs 1/L 
give 1/L ^ values that are well fitted by 



for T < Tc, and a value of Tc that is well within errors of 
(the value we find below) Tc — 0.86±0.07. In addition, 52 
and qi extrapolate to roughly the same value, for T < Tc. 
This would be in accordance with the droplet model of 
spin glasses. However, for reasons given below, this is 
not a firm conclusion. 

In principle, the critical exponent rj can be obtained 
from the plots of q2 vs 1/L shown in Fig. [31 making 
use of q2 ~ l/L'^"^"'"^ at Tc, which follows from finite 
size scaling i^i^i^ In fact, however, no meaningful num- 
ber was obtained for 77, because the errors turned out to 
be too large. 

Similar plots of q2 and qf vs 1/L for 2D are shown in 
Fig. m They clearly suggest that, at least for T > 0.4, 
92 ^ as 1/L 0. 

In order to examine the data we have for qi and q2 in 
a slightly different way, we define. 

Note that U12 is scale free, and is consequently only a 
function of ^/L, according to finite size scaling (FSS) 
theory . ^'^'^1'^^ M12 is analogous to Binder's ratio it24, which 
is defined in terms of 54 and (72 -"^^ Clearly, U12 — 1 for 
the droplet model. On the other hand 1*12 ~ for a 
macroscopic paramagnetic system, since q is normally 
distributed then, as follows from the central limit theo- 
rem and the fact that ^ is finite in a paramagnet. 

Replacement of qk by qk for all k in Eq. ([7]) gives the 
definition of Ui2. 

From TMC simulations we obtain the equilibrium re- 
sults for the RAD model in 3D that are shown in Fig. 
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FIG. 5: (Color online) ui2 vs T for systems of L x L x L spins 
in 3D. T, •, and ■ are for L — 4,6, 8, and 12, respectively. 
In addition, data points (o) for ui2 are given for L = 12. 
Lines are guides to the eye. Error bars are roughly given by 
the icon sizes. 




FIG. 6: (Color online) Plot of Ui2 vs T for 2D systems of 
L X L spins, for L = 4 (•), L = 8 (♦), and L = 16 (■). Open 
icons (o) stand for ui2 for L = 16. Lines are guides to the 
eye. Parameters for all (TMC) runs are given in table I. Icons 
are approximately of error bar size. 

[51 Data points for ui2 are also shown for L = 12 in Fig. 
[5] in order to illustrate the kind of agreement we obtain 
between ui2 and ui2. It is not clear in Fig. O whether 
Ui2 becomes approximately independent of L or keeps in- 
creasing with L for larger values of L and T < 0.9. Size 
independence then imphes 52/91 — 1 + 0.3T for T < Tc- 
On the other hand, M12 ^ 1 as L — s- 00, would give 
q2 ~ Qi for macroscopic sizes, which would be in agree- 
ment with the tentative inference we drew from Fig. [31 
We are thus led to 

1 < ^ < 1 + 0.3r (8) 

for macroscopic sizes, which does not discriminate be- 
tween the Droplet and RSB pictures of the RAD model. 

For T > 0.9 we have plotted (not shown) U12 vs 1/L, 
using data points from Fig. \5[ Such plots point to 
U12 — > as 1/L 0, which in turn implies there is a 
paramagnetic phase for T > 0.9. 



It is interesting to compare the above results with the 
ones we obtained for 2D. Plots of uu vs T are shown in 
Fig. [6l for various system sizes. Data points for U12 are 
also shown for L = 16. In contrast with the results for 
3D, the three curves in Fig. [6l appear to come together 
only gradually, as T — > 0. Plots (not shown) of U12 vs 
1/L, can be made from the the data points shown in 
Fig. [51 One can then extrapolate U12 to 1/L ^ 0. At 
least for 0.2 < T, U12 — > 0, which is consistent with a 
paramagnetic phase. 

We can obtain ^ (of an infinite size system) making 
use of the data for M12 and of the fact that, according 
to FSS^i^ ui2 is only a function of ^/L. Note that 
^/L is constant for any horizontal line that intersects the 
all the curves in either Fig. [51 or Fig. [51 We can thus 
obtain, ^(r„)/L„ = c, where c is some constant and T„ 
and L„ are the values of T and L where a horizontal 
line crosses the nth curve in Figs. [51 or [51 Different 
horizontal lines give different values of c which can be 
chosen independently in order to collapse all plots of ^ vs 
T into a single ^ vs T curve. Thus, we obtain the plots 
shown in Figs. [3 and [51 




FIG. 7: Plots of cC^^"" vs T for 3D and the values of u that 
are given in the graph, c is some undetermined constant. 
Lines are guides to the eye. 

Extrapolations in plots such as the ones shown in Fig[71 
give Tc ~ 0.88 for 3D. From the errors in the data for 
M12, we estimate an error STc — 0.05. We determine the 
exponent 1/, in £^ ^ [T — Tc)~'^ , from these plots. The 
value v ~ 1.5 gives the best straight line fit in the 0.88 < 
T < 1.2 range. On this basis we adopt the value ~ 1.5. 
Fits obtained from v values outside the 1 ^ < 2 range 
show significant curvature, whence we assign the error 
Si' — 0.5. Proceeding similarly for 2D, using plots as the 
ones shown in Fig. [3 we obtain Tc ~ 0, though a spin 
glass phase below T ~ 0.1 is conceivable, and v — 3ztl. 

From different extrapolation procedures we have ar- 
rived at values of Tc in the [0.83,0.88] range. Considering 
all the errors involved, we arrive at 

= 0.86 ±0.07 (9) 

for the RAD model in 3D. 
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FIG. 8: Plots of c^-^/" vs T for 2D the values of v that are 
given in the graph, c is some undetermined constant. Lines 
are guides to the eye. 



IV. CONCLUSIONS 



formation about critical behavior should be drawn from 
this equation, because it is not sufficiently accurate for 
it. From extrapolations of q2 and to the l/L ^ limit 
(see Fig. [3]), one might be tempted to infer that q2 and 
q\ become then equal, as in the droplet model. However, 
plots of ui2 vs T, shown in Fig. [5l do not provide firm 
support for such a conclusion, because the L oo limit 
of ui2 in the spin glass phase seems uncertain. We can 
only be reasonably sure that the limit is somewhere be- 
tween the value of Ui2 shown for L = 12 and 1. From 
this, Eq. ([5]), which does not discriminate between the 
applicability of the droplct^^ or RSB^^ pictures to the 
RAD model, follows. Results for the correlation length 
^, exhibited in Fig. [71 are consistent with ^ ~ (T — Tc)"", 
where i/ = 1.5 ± 0.5^ 

In 2D, the paramagnetic phase covers the T > range, 
though we cannot rule out a spin glass phase below T ~ 
0.1. Results for ^, exhibited in Fig. [51 are consistent with 
^ ~ T-", where u = 3±l. 



In sum, we have studied the equilibrium behavior of 
the RAD model by means of tempered Monte Carlo sim- 
ulations. The sizes of the systems we have simulated, 
temperatures, as well as other parameters, are given in 
Table I. From them, we have drawn quantitative evi- 
dence that points to the following conclusions. In 3D, 
the paramagnetic phase covers the T > Tc range, where 
Tc = 0.86 ± 0.07. For T < T^, there is a spin glass phase. 
In it, the overlap parameter, defined in Sec. Ill C[ does 
not vanish. It is approximately given by Eq. ([6]). No in- 
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